Flux Pinning and Phase Transitions in Model High- Temperature 
Superconductors with Columnar Defects 



K. H. Lee and D. Stroud 

Department of Physics, Ohio State University, Columbus, Ohio 43210, USA 

S. M. Girvin 

Department of Physics, Indiana University, Bloomington, Indiana 47405, USA 

(February 1, 2008) 

We calculate the degree of flux pinning by defects in model high-temperature superconductors 
(HTSC's). The HTSC is modeled as a three-dimensional network of resistively-shunted Josephson 
junctions in an external magnetic field, corresponding to a HTSC in the extreme Type-II limit. 
Disorder is introduced either by randomizing the coupling between grains (Model A disorder) or by 
removing grains (Model B disorder). Three types of defects are considered: point disorder, random 
line disorder, and periodic line disorder; but the emphasis is on random line disorder. Static and 
dynamic properties of the models are determined by Monte Carlo simulations and by solution of 
the analogous coupled overdamped Josephson equations in the presence of thermal noise. Random 
line defects considerably raise the superconducting transition temperature T C (B), and increase the 
apparent critical current density 3 C (B,T), in comparison to the defect-free crystal. They are more 
effective in these respects than a comparable volume density of point defects, in agreement with the 
experiments of Civale et al. Periodic line defects commensurate with the flux lattice are found to 
raise T C (B) even more than do random line defects. Random line defects are most effective when 
their density approximately equals the flux density. Near T C (B), our static and dynamic results 
appear consistent with the anisotropic Bose glass scaling hypotheses of Nelson and Vinokur, but 
with possibly different critical indices: transverse correlation length exponent u± ~ 1.2, anisotropy 
exponent z ~ 1.4 ± 0.2 (where z is defined by fii cx , with £m and £± are the correlation lengths 
parallel and perpendicular to the flux), and dynamical critical exponent z' « 2.0. 

PACS numbers: 74.50.+r, 74.60.Ge, 74.60.Jg, 74.70.Mq 
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I. INTRODUCTION. 

A major problem restricting the practical use of high- 
temperature superconductors (HTSC's) is the difficulty 
of producing a large critical current, especially in a mag- 
netic field. 1 Much of this difficulty is thought to result 
from dissipation due to flux motion - a dissipation gen- 
erally known at low or high dissipation rates as "flux 
creep" 2-6 or "flux flow" 7 respectively. When a current 
density J is introduced into the HTSC, it produces a 
J x B force (known as a Magnus force) on the flux lines. 
This force tends to set the lines in motion, producing re- 
sistive dissipation, unless appropriate defects, known as 
pinning centers, can prevent this motion, or at least raise 
the current density at which it begins. 

Recently, Civale et al 8 in an elegant set of experiments, 
have shown that columnar defects, introduced parallel to 
the flux lines by heavy ion irradiation, can greatly in- 
crease the critical current at which flux motion dissipa- 
tion begins, relative to the point defects which arc more 
commonly introduced as pinning centers, e. g. by pro- 
ton irradiation 9-11 The same columnar defects were also 
found to increase the temperature of the so-called "irre- 
versibility line" 12 in the magnetic field-temperature (H- 
T) plane, below which flux motion essentially ceases in 
the limit of a weak applied current. The columnar pins 
were produced by irradiating the HTSC with a beam of 
heavy ions parallel to the c-axis. It is not surprising that 
columnar defects should be effective pins: they provide a 
long pinning center which should provide a much stronger 
pinning potential for a long flux line than will an equal 
concentration of point defects. However, a realistic cal- 
culation which demonstrates this effect has been lacking. 

In this paper, we present some simple model calcula- 
tions which demonstrate both of the effects observed by 
Civale et al 8 and also suggest some alternative meth- 
ods for further increasing both the critical current and 
irreversibility temperature of HTSC's. Our approach 
is to describe the HTSC as a three-dimensional collec- 
tion of resistively-shunted Josephson junctions (RSJ's), 
in which temperature is simulated by a Langevin noise 
source of the appropriate strength in each junction. 13 
Such a model is obviously far from a realistic HTSC. 
However, the model does contain some of the essential 
physics needed to describe transport in HTSC's: it em- 
bodies coupled fluctuating phases within the context of 
a reasonable dynamics, and it allows for the introduc- 
tion of an applied magnetic field in a simple way. In 
this view, the Josephson-coupled "grains" should proba- 
bly be considered as representing small patches of phase- 
coherent superconductor, of dimensions comparable to 
the coherence length. 14 Thus, the model is not restricted 
to literally granular materials, but could apply to single 
crystals with more microscopic disorder, in the extreme 
type-II limit (penetration depth A much larger than co- 
herence length £).. We have shown elsewhere that a simi- 
lar model describes the difference between transverse and 



longitudinal magnetoresistance of a HTSC, in qualitative 
agreement with experiment. 15 

We turn now to the body of the paper. Section II de- 
scribes the model for both the thermodynamic and trans- 
port properties of the HTSC with defects. Section III de- 
scribes our numerical results for these properties. A brief 
discussion follows in Section IV. Three Appendices sum- 
marize the static and dynamic scaling hypotheses used 
to analyze our numerical results. 

II. MODEL. 

A. Thermodynamics. 

We consider a simple cubic three-dimensional network 
of N superconducting "grains" weakly coupled together 
by Josephson junctions, and driven by an externally ap- 
plied current. The \ th "grain" is described by a supercon- 
ducting order parameter \ipi\ exp(i0j). We neglect fluctu- 
ations in the amplitude \ipi\ and allow the phase Oi to 
fluctuate. With these assumptions, the thermodynamic 
properties of the network are given by the Hamiltonian 

H = - J2 E J;a cos (^ - 6 J - A a) (!) 

<ij> 

where 

= ~2^ c > i 3 (^) 

is the coupling energy between grains i and j, l C ;ij is 

the 

corresponding critical current, 

Aii= $~oL A ' d1, (3) 

$o = hc/2e, and A is the vector potential, taken to 
be that of the externally applied magnetic field (this is 
equivalent to assuming a Josephson penetration depth 
large compared to the intergranular separation). The 
sum runs over distinct nearest neighbor pairs. 

Given H, equilibrium properties are obtained via an 
average with respect to to a canonical ensemble. Thus, 
for example, the average of some operator O(0i, ...On) is 
obtained from 

<o>= J O(0 u ...,o N ) e - H ^--' e ^/ kBT n^ 1 d0 t /z 

(4) 

where 

Z = J UidOiexp(-H/k B T). (5) 
is the canonical partition function. 
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In the calculations to be described below, we have gen- 
erally dealt with disordered samples. In that case, we cal- 
culate averages both over a canonical ensemble (denoted 
< ... >) and over different realizations of the disorder 
(denoted [...]). We have considered primarily the specific 
heat per grain Cy and the so-called hclicity modulus ten- 
sor with components 7y . Cy is generally computed from 
the fluctuation expression 



Cy = [< H 2 > - < H > 2 }/{Nk B T 2 ). 



(6) 



The helicity modulus (or equivalently, the supcrfluid 
density) is the free energy cost of imposing a twist in 
the phase at the boundaries of the sample. Its princi- 
pal elements are essentially spin-wave stiffness constants. 
Rather than imposing a twisted boundary condition and 
calculating the resulting increase in free energy, it is more 
convenient to use periodic boundary conditions and cal- 
culate 7^ as 



7»j 



o 2 f 
a Aid a'. 



(7) 



A'=0 



Here A' represents an added uniform vector potential (in 
addition to that which produces the applied magnetic 
field) which is included in the Hamiltonian in order to 
produce a twist. The various second derivatives in eq. 
(7) are readily computed explicitly for an ordered or a 
disordered sample, with the result for, e. g., 



k B T 



( Y. Ej-sitfj coB(6i - 6j - Ay) 

\<ij> 

^ Ej-ijXij sin(0i - 9j - Ai 



><y> 



Xh sin 



9j Aij ) 



\<ij> 



(8) 



where the x coordinate of the distance 

between nearest neighbor grains i and j . Similar expres- 
sions hold for the other components of -y. 16 ' 17 



B. Dynamics. 

There are many dynamical models whose equilib- 
rium thermodynamic properties are represented by the 
model just described. We choose a dynamical model 
corresponding to a simple cubic array of overdampcd 
resistivcly-shunted Josephson junctions, driven by an ap- 
plied current. The network is then characterized by the 
set of coupled equations 



Iij = I c sin(0i - 6j -Aij) + ^-+ I L . t ij 
V^V t -V^^±(0 t -0 3 ) 



E 7 > 



ij 1 i:ext 



A-dl. 



(9) 

(10) 
(11) 

(12) 



Here is the current from grain i to grain j, which is 
written as the sum of a Josephson current and an Ohmic 
current through the shunt resistance Ry ; Vy is the volt- 
age difference between grains i and j; and \i- ex t is the 
external current fed into grain i. In the calculations de- 
scribed below, the current is always fed into one face of 
the array (an equal amount I into each grain) and ex- 
tracted from the opposite face, with periodic transverse 
boundary conditions. Eq. (11) is Kirchhoff's law describ- 
ing current conservation at grain i. lL-ij is a Langevin 
noise current 18 between grains i and j, introduced to sim- 
ulate the effects of temperature, which satisfies the rela- 
tion 



<lL;ij(t) >=0 



< lL;ij(t)lL-Mt') > = 



2k B T 
FL 



Sij;klS(t - 1'), 



(13) 
(14) 



where T is the absolute temperature and the brackets 
denote an ensemble average. We solve these equations 
by Euler iteration, as described previously. 13 ' 15 



C. Geometrical Models for Line and Point Disorder. 

We have considered two types of models to describe 
disorder, which we denote Model A and Model B. In 
Model A, the bond energy Ej-ij between grains i and 
j is assumed to vary randomly between and twice its 
average value. In Model B, we introduce disorder sim- 
ply by removing a certain fraction of the grains, as well 
as their associated Josephson junctions (but not shunt 
resistances). 

For either Model A and Model B, we can assume ei- 
ther "line disorder" or "point disorder." In Model A, 
point disorder can be introduced by assuming that the 
bond strengths of different bonds are completely uncorre- 
cted. "Line disorder" consists of assuming that the bond 
strengths are uncorrelated in the xy plane, but are per- 
fectly correlated in the z direction - that is, the strength 
of a given bond, whatever its orientation, depends the x 
and y coordinates describing its location, but not on the 
z coordinate. 

To introduce point disorder in Model B, we remove the 
grains at random. For line disorder, we remove lines of 
grains parallel to the z axis. The removal of these grains 
effectively converts the neighborhood of the grain from 
superconducting to normal. In Model B, we have also 
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considered "periodic line disorder," in which line defects 
are arranged periodically in the xy plane, as described 
further below. 



D. Magnetic Field. 

We consider magnetic fields applied in the z-dircction, 
i.e., parallel to the line defects. In previous calculations of 
this sort, it has been standard to use the Landau gauge, 
A = Bxy. This gauge severely restricts the possible mag- 
netic fields that can be considered if one also requires pe- 
riodic boundary conditions in all three directions (as in 
the Monte Carlo simulations) or in the transverse direc- 
tions (as with the dynamic simulations). We therefore 
use a different gauge previously used by Arovas and Hal- 
danc in other contexts 19 . 

To define this gauge, we consider an L x L square array 
of lattice constant a, with the origin taken as the lower 
left hand corner of the array. Then we take Ay = 2wfn 
for bonds pointing in the y direction and located at x = 
na; Ajj = for all bonds in the x direction except for 
those in the extreme right-hand column of plaqucttcs, 
and Aij — -2wfLm for horizontal bonds in that extreme 
right-hand column, at y = ma, x=La. (In this expres- 
sion f = &/$o, where $ is the flux per plaquette.) With 
this choice of gauge, it is readily verified that the fac- 
tors sum to 2nf (modulo 2tt) around each plaquette, 
as required. The requirement of periodicity in the two 
transverse directions will now be satisfied as long as f is 
a multiple of 1/L 2 . This is much weaker than the con- 
dition imposed by periodicity when the Landau gauge is 
used, which is that f be a multiple of 1/L. Hence, a much 
finer grid of magnetic fields can be considered than in 
most previous calculations. 



III. NUMERICAL RESULTS 
A. Model A. 

We begin by presenting our numerical results for Model 
A. Fig. 1 shows the specific heat Cy per grain and the 
helicity moduli 7|| and 7j_ parallel and perpendicular to 
the magnetic field, for an ordered array of size L x L x L 
grains with periodic boundary conditions, isotropic cou- 
pling, magnetic field f = 1/4, and several different sizes 
(L = 8, 12, and 16). Cy shows clear signs of diverging 
in the vicinity of k^T = l.lEj (the peak in Cy is grow- 
ing with increasing lattice size, suggesting a continuous 
phase transition). We interpret this temperature as a 
melting transition with a discrete symmetry associated 
with the periodic lattice. Similar numerical results have 
been obtained for this model by Shih et al 17 ' 20 . 

We may get a clearer quantitative picture of this phase 
transition by applying a static scaling analysis, 21 ~ 23 ' 16 as 
described in detail in Appendix A. In Fig. 2, we plot L711 



and L7^ for several values of L as a function of tempera- 
ture. As is clear from the Figure, each attains a universal 
value near the same temperature k^T c = l.lEj, suggest- 
ing that this is indeed the critical temperature for this 
model. When combined with the scaling analysis of Ap- 
pendix A, these results suggest that this phase transition 
is characterized by a single correlation length £, i.e. z=l, 
where z is the anisotropy exponent, defined by the rela- 
tion £|| cx where £|| and £j_ are the correlation lengths 
parallel and perpendicular to the field. 

Fig. 3 shows the quantities of Fig. 1, but with line 
disorder, 24 Once again, we consider isotropic coupling 
strengths and use f = 1/4, and several different (but cu- 
bic) box sizes. The calculations shown are the result of 
averages over many realizations of the disorder, as indi- 
cated in the legends of the Figure. In contrast to Fig. 
1, there is very little size dependence of Cy, suggesting 
that Cy either does not diverge or at most diverges very 
weakly. Note also that the helicity modulus 7|| goes to 
zero at a substantially higher temperature than in the 
ordered case, indicating that the superconducting tran- 
sition temperature is increased by the line disorder. 

Another point is that, like the ordered case, 711 also 
seems to vanish continuously with temperature. If we 
assume a power law of the form 7|| cx \T — T c | 9 " , then 
Fig. 3 suggests g|| « 0.5 and k B T c « 1.7 < Ej > The 
the numerical uncertainties in -f± are much larger, so an 
accurate estimate of the analogous quantity g± is diffi- 
cult. Indeed, at any given temperature, the Monte Carlo 
convergence of j± is much slower than for 7||. This is 
presumably because our model is both frustrated and 
disordered in the xy plane, but is neither frustrated nor 
disordered in the z direction. Hence, the system rapidly 
responds to any twist in that direction, but much more 
slowly in the xy plane. Also, j± converges more slowly 
in the disordered case than in the frustrated but ordered 
model of Fig. 1, suggesting that the slow convergence 
is caused by the huge number of metastable states of 
nearly equal energy which are expected in the disordered 
system. 

From the static scaling analysis of Appendix A, we 
can roughly estimate the anisotropy exponent z defined 
there. 21 From cq. (23) of Appendix A, = v±(2 — z), 
where u_\_ is related to the transverse correlation length 
£± by £_l cx |(T - T c )/T c \- V± . Hence, our numerical 
results suggest 

i/i(2-z)rj0.5. (15) 

If v±_ is finite, this equality suggests that z < 2 for this 
model. Secondly, since Cy is apparently nondivergent, 
cq. (21) shows that 

{2 + z)v^>2. (16) 

Adding these results, we get 4z/j_ > 2.5 or v± > 0.63. 
Hence, (2 - z) = 0.5/v± < 0.8 or z > 1.2. Combining all 
these arguments, we estimate 1.2 < z < 2. 
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In order to go further, we must estimate g±. According 
to eq. (24) of Appendix A, g± = zv^. From Fig. 3, we 
estimate zz/j_ w 1 — 2. Combining this with eq. (15), 
and assuming zz/j_ = 1.5, we get z « 1.5, v± w 1, which 
satisfy the inequalities described above. They also agree 
with Chayes et aP 5 , who have proposed for a wide class 
of continuous-spin models with disorder, that v± > 1 
rigorously. Our estimate is, of course, subject to large 
numerical uncertainties, especially in gj_. 

We turn now to the dynamical properties of Model A, 
concentrating on isotropic coupling with and without line 
disorder. Fig. 4 shows the IV characteristics of a 6 x 6 
x 9 array at magnetic field f = 1/4, with no defects and 
with current density J _L B, plotted at several different 
temperatures. Figs. 5(a) and 5(b) show the IV charac- 
teristics for the analogous model with line disorder, a 6 x 
6x6 unit cell, and two current orientations: J _L B and 
J\\B. Like the Monte Carlo results, the IV characteris- 
tics also suggest that T c is increased by line defects. To 
make this clearer, we have plotted in Fig. 6 the resistivity 
p =< V > /(LRI) at a current level I=0.05I C : the in- 
troduction of line defects reduces p at every temperature, 
relative to the no-defect case. 

The shape of the IV characteristics is also changed by 
the introduction of line defects. With no line defects and 
J _L B (Fig. 4), there is a fairly clear critical current 
onset for temperatures fc^T < LOEj/ks- When line 
defects are present, the IV characteristics suggest no clear 
critical current for J _L B (the analytic form of these IV 
characteristics is discussed further below). By contrast, 
for J\\B, a critical current seems to develop for UbT < 
1.3 - lAEj. 

We have attempted to scale the IV characteristics of 
Fig. 5 according to the formalism 21 outlined in Ap- 
pendix B. For J _L B, we plot F,±t~ v± ( 1+Z ' against 
J±t -v±(i+z) ( w here t = |T - T c \/T c ) for various esti- 
mates of T c , z, z',and i> ± . For J\\B, we plot F,\\t- V±( - Z+Z "> 
against J\\t~ 2v± . Our best results are shown in Figs. 
7(a) and 7(b). For a given choice of T c , the best fits 
seem to correspond to somewhat different values of the 
parameters in the parallel and perpendicular directions 
- not surprising in view of the numerical uncertainties, 
small sample sizes, and limited current ranges. For 
J _L B, assuming ksT c = 1.7Ej, we get v±(l + z) w 3.8, 
+ z') rts 4.5. For J\\B, on the other hand, we get 
2v± « 2.4, v±(z + z') « 3.4. Combining these results 
with the previous Monte Carlo estimates, we estimate z 
= 1.4 ±0.2, z' = 1.8 ±0.3, and v ± = 1.3 ±0.3. The error 
estimates in all three cases are simply subjective assess- 
ments of our confidence in these numbers over the ex- 
tremely limited current and size ranges considered. This 
estimate is also based on the assumption that the same 
value of z' applies to both the parallel and perpendicu- 
lar case. That assumption is apparently not valid in the 
case of short-range interactions between the vortex lines, 
as has recently been shown by Wallin and Girvin. 26 



Both above and below the assumed T c , the IV charac- 
teristics in both the parallel and perpendicular case col- 
lapse adequately (though not perfectly) onto universal 
scaling functions above and below T c . For both current 
directions, but especially for J\\B, there is a conspicuous 
ohmic tail in the IV characteristics below T c . We believe 
this is a finite-size effect, as further analyzed in Appendix 
C. Indeed, we have checked that the tail becomes smaller 
and smaller as the size of the array is increased. When 
T < T c , the perpendicular and parallel scaling functions 
are quite dissimilar. In both cases, we can obtain ac- 
ceptable fits with scaling functions of the form E(x) oc 
exp[—(A/x)^], but the best fit for fi± is in the range 
0.2-0.4 while fj,^ »1. This is consistent with the IV char- 
acteristics of Fig. 5, which show a much clearer critical 
current developing in the parallel direction than in the 
perpendicular direction. 

B. Model B Disorder. 

To study Model B disorder, we consider an 8 x 8 x 5 
lattice with flux per plaquettc f = 1/8, parallel to the z 
(thin) direction. We have considered three types of defect 
configurations: (i) "point defects," introduced by remov- 
ing at random 40 grains, and their associated Josephson 
couplings (but not shunt resistances) from among the 320 
grains of the lattice; (ii) "random line defects," consisting 
of eight line defects, each five grains long, parallel to the 
z direction but randomly distributed in the xy plane; and 

(iii) "periodic line defects," in which the line defects are 
arranged with the periodicity of the ground state phase 
configuration at f = 1/8. For reference, we also consider 

(iv) the no-defect configuration. 

In the absence of an applied current, the phases of 
configuration (iv) will settle into a z-independent ground 
state configuration. This can be found numerically, e. 
g., by starting the phases in a random arrangement and 
iterating the Josephson equations at zero applied current 
until a state of no voltage is obtained (care must be ex- 
ercised to avoid falling into a metastable minimum). To 
calculate the IV characteristics, we typically begin with 
this ground state, gradually increasing the applied cur- 
rent for various defect configurations. 

Fig. 8 shows the resulting IV characteristics (for 
J ± B) at temperature T = 0. Case (iv) has a crit- 
ical current w 0.12I C per junction, comparable to the 
calculated depinning critical current for a single vortex 
in a large square array 27 . This suggests that the crit- 
ical current is not too much influenced, in this case, 
by vortex- vortex interactions. The critical current is 
increased slightly by point defects [case (i)], somewhat 
more by random line defects [case (ii)], and considerably 
more again by periodic line defects [case (iii)] . In case (ii), 
the functional form of the IV characteristic is consider- 
ably modified by the defects (being concave up rather 
than concave down). We find that it is fairly well fitted 
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by < V > /{LRI C ) = Aexp[-C(I c /iy} with /x » 0.3, A 
« 3 x 10 4 , and C « 10. 

Fig. 9 shows the temperature-dependent resistivity 
p(T) = V/I at a current of 0.1I C per junction, for cases 
(i) - (iv) and J _L B. For reference, the superconduct- 
ing transition temperature T C (B = 0) of the Hamilto- 
nian (1) in zero magnetic field and zero current is known 
to occur at k B T c « 2.21_E,/. 28 Hence, Fig. 9 suggests 
that, whatever the defects, T c (f = 1/8) < T c (/ = 0) - 
that is, as expected, T c is reduced in a magnetic field. 
However, relative to the zero-defect case, T c (l/8) is in- 
creased slightly by point defects, more yet by random line 
defects, and still more by periodically arranged line de- 
fects. The increase in T c (l/8) produced by random line 
defects corresponds to the increase in the "irreversibil- 
ity temperature" observed by Civale et al when random 
line defects are introduced parallel to the magnetic field. 
Note that, in Fig. 9, the density of line defects "equals" 
the density of point defects in the sense that an equal 
amount of superconducting material is removed in each 
case. The only difference among the various defect curves 
in Fig. 9 is the degree to which the disorder is correlated. 
Hence, this plot provides a very direct illustration of the 
influence of correlation in raising the "irreversibility tem- 
perature." 

Fig. 10 shows the T = IV characteristics for a flux 
density f = 1/8 and several densities f^ of randomly dis- 
tributed line defects. Fig. 11 shows a similar plot for p(T) 
at a current level of 0.1 I c . These Figures suggest that 
(at least for f = 1/8) both the T = critical current and 
T C (B) are largest when f « f^, a conclusions which may 
qualitatively agree with experiment 8 . To confirm these 
conclusions, however, Monte Carlo calculations should be 
carried out on the analogous Hamiltonian to determine 
the dependence of T c on defect line density. 

To summarize for Model B disorder, our results show 
that columnar defects oriented parallel to the flux lines 
tend to increase the critical current, and to push up the 
superconducting transition temperature T C (B), relative 
to the same number of random point defects at the same 
field, in apparent agreement with experiment. We pre- 
dict also that a periodic arrangement of line defects com- 
mensurate with the defect-free flux lattice will be even 
more effective in increasing both the critical current and 
T C (B). Finally, we have preliminary evidence that the 
pinning (and T c -enhancing) effects of random line defects 
are optimized when the defect density is comparable to 
the flux density. 

IV. DISCUSSION. 

Nelson and Vinokur 21 have recently proposed a the- 
ory of the superconducting transition in materials with 
correlated disorder. When the density of line defects is 
greater than the density of flux lines, their theory leads 
to a phase transition from a high-temperature flux liq- 



uid into low-temperature "Bose glass" phase. In their 
theory, the three-dimensional superconductor maps onto 
a two-dimensional system in which the flux lines play 
the lines of interacting Bose particles while the line de- 
fects become static point defects. When the density of 
line defects equals the density of flux lines, Nelson and 
Vinokur predict instead a transition into a "Mott insu- 
lator" phase, in which the flux lines are localized on the 
line defects. 

Although our model is presumably not in the same 
universality class (because of the long range interactions 
among the vortices) as that considered by Nelson and 
Vinokur, we briefly interpret our results in the context of 
this theory. The "Bose glass" regime, in which the line 
defect density id > /, corresponds to our Model A and 
one of the cases considered in our Model B. In this regime, 
for T < T C (B), Vinokur and Nelson predict a nonlinear 
IV relation of the form V/J cx cxp[-(E k /k B T)( J / J) 1 / 3 ] 
where J is the current density and Ek and J are con- 
stants depending on the strength of the pinning poten- 
tial and other parameters. Our IV characteristics for 
Model A (with JIB) are not inconsistent with this 
behavior, though our samples are extremely small in 
the dynamical calculations. (Although the number of 
flux lines is very small in these calculations, the num- 
ber of flux loops, formed by localized bowing out of the 
flux lines in response to applied currents, is consider- 
ably larger. Since dissipation occurs largely by motion 
of these loops, in the Bose glass model, our dynamical 
samples may not be as inadequately small as appears on 
first sight.) For J _L B, a simple activation form [V/J 
oc exp(— Ek Jo/( J&bT))] seems ruled out by our results. 
For J\\B, our data seems consistent with an activated 
behavior of this form; however, we are not aware of an 
analytical theory describing Bose glass transport in this 
regime. Obviously, more detailed numerical simulations, 
involving much larger numbers of line defects and more 
disorder realizations, are needed before any definite con- 
clusions can be drawn about the quantitative applicabil- 
ity of the Bose glass picture to this dynamical model. 

Our results for Model A disorder seem consistent with 
both static and dynamic scaling hypotheses as applied 
to this anisotropic phase transition, with critical indices 
z«1.4± 0.2, z' w 1.8, v±_ w 1.3 ± 0.3. z is, as expected, 
smaller than that found in analogous calculations with 
short-range interactions 26 between vortex lines, for which 
z = 2, but is, perhaps coincidentally, in the range of find- 
ings for other lattice models with line disorder 29 . Our re- 
sults (although based on very small dynamical samples) 
seem to suggest there may be only a single dynamical 
exponent z' for transport both perpendicular and par- 
allel to the line defects, in contrast to the short-range 
model. This conclusion may, however, be a function of 
the particular dynamics assumed. 30 The value of v±_ is 
comparable to values of v which have been reported in 
the experimental literature 31 in cases where point disor- 
der is expected to dominate; to our knowledge, no direct 
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measurement of v±_ has been carried out for samples with 
known line disorder. 

To summarize, we have studied flux pinning by defects 
in three-dimensional Josephson-junction networks with 
various types of point and line defects. We find that 
line defects considerably enhance both the critical cur- 
rent and the upper critical field, relative to the same con- 
centration of point defects. We also find some indications 
that the normal-to-superconducting transition in the case 
of random line disorder is a continuous phase transition 
marked by both static and dynamical critical phenomena. 
These conclusions are, however, subject to large numeri- 
cal uncertainties arising from the relatively small samples 
studied. The resulting critical exponents may be appli- 
cable to high-T c superconductors in which the intrinsic 
anisotropy is not too great (such as YBa 2 Cu307_5) and if 
screening fields can be neglected (extreme Type II limit). 
Larger scale calculations, including screening, and mak- 
ing use of more realistic pinning models, will be neces- 
sary for quantitative comparisons to experiment. These 
are planned for future work. 
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Appendix A. Static Scaling. 

We describe our numerical results for line defects 
within the framework of a scaling analysis suitable for 
both our static and dynamic results, based largely on 
previous discussions of Nelson and Vinokur. 21 In this Ap- 
pendix, we describe the static scaling hypotheses. 

Consider a HTSC with line defects and a magnetic field 
both oriented in the z direction. Suppose that there is 
a phase transition at some temperature T C (B), where B 
is the magnitude of the applied magnetic field. Assume 
also that this phase transition is characterized by two di- 
verging correlation lengths, £j_ and £||, corresponding to 
correlations in the xy plane and z direction, respectively. 
To allow for the possibility that these diverge with dif- 
ferent critical exponents, we write 



f ll oc r ZUl - 

t=\T-T c (B)\/T c (B). 



(17) 
(18) 
(19) 



An isotropic phase transition is a special case of this be- 
havior with z = 1. 

By analogy with the usual isotropic hypcrscaling ex- 
pression for the singular part f s of the free energy density 
near T c , we assume that i s behaves as 



Pis OC 



^11 



(20) 



(where /3 = l/fe^T). Hence, the specific heat has a sin- 
gularity of the form 



(21) 



To estimate the behavior of 7_i_ and 7|| , the principal 
components of the helicity modulus tensor, we extend an 
argument of Cha et aP 3 to anisotropic phase transitions. 
First imagine that the array is subjected to a phase gra- 
dient \7 Z 8 in the z direction. The change in free energy 
per unit volume is 



/^oc^l^ocj, 



(22) 



where we have replaced V z by the inverse of the charac- 
teristic length With the use of eqs.(18) and (20), this 
gives 



7|| o, jj. K t <a-.W 



(23) 



where the last proportionality describes the expected 
critical behavior near T c . A similar argument applied 
t° 1± gives 



7_L oc \- oc t zu± - . 

qi 



(24) 



In a Monte Carlo calculation, it is necessary to calcu- 
late these quantities in finite-size sample, usually a par- 
allelopiped of volume Lj_L||. The natural scaling form 
for the helicity moduli in such samples is 



7|| = #G^||/L||,a/i±), 



(25) 
(26) 



where F(u, v) and G(u, v) are universal functions. Ex- 
pressions (25) and (26) can be written in more conve- 
nient forms by making the change of variables F(u, v) = 
uH(uv~ z ,u), G(u,v) = (v 2 /u)K(uv~ z ,u) to yield 



^ = —H{L z jL^m^) 



in 



T|| = jr K ( L 'd L h L \\IZ\\)- 



(27) 
(28) 



At T = T c , the second argument of both H and K van- 
ishes. Hence, for a finite sample whose dimensions are 
chosen such that Lm oc L z ,, it follows that 



£||7_L = const.; 
Lj_7||/L|l = const. 



(29) 
(30) 
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These relations can, in principle, be used to determine 
the transition temperature with high accuracy by exam- 
ining the behavior of the components of 7 in a series of 
boxes of different volumes, such that the ratio Lm/Lj_ is 
held constant. The method is to plot L||7j_ and Lj_7||/_L|| 
for different volumes; all should cross at T = T c . Un- 
fortunately, this method works only provided z is known. 
Since z is apparently in the range 1.2 - 1.5 for the present 
model, but difficult to determine with greater accuracy, 
we have not attempted this kind of anisotropic finite-size 
scaling in the present paper. 

Appendix B. Dynamic Scaling. 

For dynamical quantities, we may again follow 
and somewhat extend the arguments of Nelson and 
Vinokur. 21 We consider first the electric field Ej_ and 
current density Jj_ in the transverse direction. In this 
case, we postulate a scaling relation of the form 

E± = C 1 a E ± ^{T l J^^/{2ek B T)) (31) 

where E± j_ are scaling functions which apply respec- 
tively above and below T c . To determine a, suppose 
the HTSC is in the Ohmic regime at T > T c . In this 
regime, E +j j_(a;) cx x, whence cr± = Jj_/E± cx £" -1-z . 
However, we also expect that o~± should scale like 7_l/w 
(that is, both of these quantities should have the same 
power-law dependence on where uo is a characteris- 
tic frequency). In turn, we expect that u> should vanish 
near this continuous phase transition, with a character- 
istic temperature-dependence given by where z' is 
a new dynamical critical exponent. Combining all these 
relations, and using eq. (24), we obtain 

a=l + z'. (32) 

Similarly, for transport parallel to the z axis, we expect 
the scaling relation 

E ll =^ 1 »E ±Al (hJ li ej(2ek B T)), (33) 

and making arguments analogous to the perpendicular 
case, wc find 

b = z + z'. (34) 

Precisely at T = T c we expect both Ej_ and Ey to 
vary as power laws in and J|| respectively. This be- 
havior implies that at T c the scaling functions F l ±,±(x) 
and E ± n(a;) should take the forms 

E± t± (x) oc x c , (35) 
E ±tll (x) oc x d . (36) 

The exponents c and d can be determined by observing 
that £j_ and £y are infinite at T c . In order for eqs. (31) 
and (33) still to be satisfied at T c , the left and right hand 
sides must should involve equal powers of t. This leads 
to the results 



c={l + z')/{l + z) (37) 
d=(z + z')/2. (38) 

Thus, calculating or measuring the current-voltage char- 
acteristics precisely at T=T C should yield power laws 
whose slopes determine the exponents z and z'. 

The above arguments assume that there is a single dy- 
namical critical exponent z'. If instead, there are two 
such exponents z' ± and z|| describing the divergent relax- 
ation times in the perpendicular and parallel directions, 
then eqs. (37) and (38) are replaced by the relations c = 
(l+z' ± )/(l+z); d = (z+z(,)/2. 

Appendix C. Dynamic Scaling in Finite Size Sys- 
tems. 

Our finite-size IV characteristics often show an Ohmic 
tail at low currents, even at temperatures well below the 
putative superconducting transition. In this Appendix, 
we give an argument suggesting that this tail is a finite- 
size effect which would disappear in sufficiently large sys- 
tems. 

We present the argument for the case J||B, where the 
numerical results most clearly show the finite-size tail. 
However, a similar argument should also hold for the 
perpendicular case. In the parallel case, for T< T c in 
an infinite system, as shown in Appendix B, the electric 
field and current density are related by 

E\\ = G iZ+Z,) E-,\\ (hJ\\e ± /(2ek B T)), (39) 

where E_ is some universal function. Now write 
E_ ||(a;) = xF_^(x). It follows that the resistivity 
Pll = E\\/ J\\ can be written as 

PI, cxC 2 -^ z >-,||(?iJ||Ci/(2efc B T)). (40) 

For a cubic system of edge L, this relation must involve 
another variable, the ratio £±/L: 

PI, ocer'~''F-,\\(hJ\\e±/(2ek B T),U/L). (41) 

The numerical data presented in Fig. 7(b) suggest that 
F_ || (x, 0) falls rapidly to zero with decreasing x. Indeed, 
the data seems to fit roughly to the relation 

F_ t \\(x,0) = F o exp(-A/ X W) (42) 

where F a and A are constants and /X|| w 1). We have 
no theory for the finite-size version of this function, 
but a plausible guess suggests itself. The dimension- 
less argument x — TiJ\\£,\ /(2ek B T) may be expressed as 
x = (£_l/£j) 2 where O = {2ek B T / (fiJ^)) 1 / 2 is a char- 
acteristic length defined by the current density Jm (note 
that length is measured in units of the intergranular sep- 
aration). When this length becomes larger than the sys- 
tem size, the current length should be replaced by L. In 
order to include this behavior in the scaling form, we may 
postulate 
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= F cxp(-A/[hJ li ej(2ek B T) + Cl/L 2 }^ )• (43) 

This function has several desirable properties. First, for 
large L, it reduces to the infinite-size form of eq. (42). 
Secondly, for very small Jm, pii becomes independent of 
J|| (i. e., becomes Ohmic) and given by 

P\\ K ^\~ Z ~ Z ' exp(-AL 2 ^" ). (44) 

Eq. (44) is the desired low-current Ohmic tail seen 
in our calculations. As expected, it goes away at large 
enough sizes, or low enough temperatures (£|| becomes 
smaller and smaller as the temperature is decreased below 
T c ). Since z + z' w 4 — 5, the prefactor in eq. (41) grows 
with decreasing temperature. However, its growth should 
be more than offset by the decreasing exponential, so that 
/dii should decrease with decreasing temperature for fixed 
large L. Thus, the argument presented in this Appendix 
gives a plausible explanation of the finite-size numerical 
results discussed in the text. 
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FIG. 1. Specific heat per grain, Cv; and parallel and 
perpendicular components of the helicity moduli, 7y and 7±, 
for an ordered L x L x L lattice of grains, with isotropic 
coupling (Ej,± = £^11 = Ej), magnetic field f = 1/4, and L 
= 8, 12, and 16, with periodic boundary conditions, plotted 
as a function of temperature T. 



FIG. 2. L7|| and L7X versus temperature, for the array of 
Fig. 1. 



FIG. 3. Same as Fig. 1, but for an array with Model A 
line disorder in the bond strengths. The calculations involve 
averages over 35, 21, and 19 realizations of the disorder for L 
= 8, 12, and 16 respectively. 



FIG. 4. IV characteristics for the model of Fig. 1, L x 
L x L z array, current density J _L B, at several different 
temperatures, averaged over three different (random) choices 
of initial conditions, and L = 6, L z = 9. In this and subsequent 
Figures, the lines are simply interpolations between calculated 
points. 



FIG. 5. IV characteristics for the model of Fig. 3, L x L 
x L array, with L = 6: (a) J _L B, and (b) J\\B, at several 
different temperatures, as indicated. 



FIG. 6. Resistivity p =< V > / (LRI) for an L x L x L 
array with L = 6, at a current level I = 0.05I C per grain, for f 
= 1/4 and no defects, J _L B (triangles); line defects, J _L B 
(squares); and line defects, J\\B (circles). Cases (b) and (c) 
each involve averages over ten realizations of the disorder. 



FIG. 7. Scaling plots of the IV characteristics from Figs. 
5 for (a) J _L B and (b) J\\B. In both cases, the IV char- 
acteristics both above and below T C (B) collapse reasonably 
well onto universal scaling functions over the limited current 
ranges considered. For T < T c , there are low-current Ohmic 
tails in both cases (especially for J\\B), which are probably 
due to finite size effects. The fitting parameters for curves (a) 
are z =1.5, z' = 2.0, v± = 1.5, ksT c = \.1Ej\ for curves (b), 
they are z=1.5, z' = 1.3, u ± = 1.2, and k B T c =1.7Ej. 
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FIG. 8. IV characteristics for Model B disorder and 
J _L B, at temperature T = 0, an L x L x L z array, with L 
= 8, L z = 5, magnetic flux f = 1/8 with B parallel to the z 
(thin)dircction: no defects (dotted curve); 40 randomly dis- 
tributed point defects (full curve, average of 7 realizations); 
8 randomly distributed line defects parallel to the z direction 
(dashed curve, average of 10 realizations); and 8 periodically 
distributed line defects in the z direction (dot-dashed curve). 
I is the applied current per grain; < V > is the time-averaged 
voltage across the sample, averaged over the directions per- 
pendicular to the current; R is the shunt resistance, and I c is 
the critical current of each junction. For reference, the critical 
current for the ordered lattice at f = is I/I c =1.0. 

FIG. 9. Resistivity p(T) = V/I, at an applied current 
I = 0.11c per junction, plotted versus temperature T, at f 
= 1/8. Symbols as in Fig. 10. Inset: ground state vortex 
line configuration for f=l/8 lattice. Filled squares denote loci 
of vortex lines (plaquettes of positive vorticity, i.e. current 
circulating counterclockwise); empty squares are plaquettes 
of negative vorticity. For reference, T c (/ = 0) m 2.2\Ej /Ub- 

FIG. 10. Same as Fig. 8, but for a flux f = 1/8 and several 
densities id of columnar defects oriented parallel to the z axis. 
Each curve represents an average over ten realizations of the 
disorder. In this case, the lattice is 8x 8 x 8. 

FIG. 11. Same as Fig. 9, but for a flux f = 1/8 and several 
densities id of columnar defects oriented parallel to the z axis. 
Each curve represents an average over ten realizations of the 
disorder. 
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